function [ALPHA,SIGMA,resids] = var_gibbs(factors,exogenous,prior,pvar,Q,DiagOpt)
%VAR_GIBBS Produce a Random Draw from the Posterior Distribution
%   of the parameters of a Vector Autoregression Model.
%   factors : vector of endogenous variables
%   exogenous: vector of exogenous regressors#
%   prior : prior specification (0 = noninformative; 1 = N-IW)
%   qq 

%--------------------------DATA HANDLING-----------------------------------

Yraw = factors;

% Get initial dimensions of dependent variable
[Traw M] = size(Yraw);

Y1 = Yraw;
Y2 = Yraw;

% Generate lagged Y matrix. This will be part of the X matrix
Ylag = mlag2(Y2,pvar); % Y is [T x M]. ylag is [T x (Mp)]
if isempty(exogenous)
X1 = Ylag(pvar+1:Traw,:);
else
X1 = [exogenous(pvar+1:Traw,:) Ylag(pvar+1:Traw,:)];
end

if size(exogenous,2)==0; constant=0; else constant = 1; end
% Get size of final matrix X
[Traw3 K] = size(X1);

% Create the block diagonal matrix Z
Z1 = kron(eye(M),X1);

% Form Y matrix accordingly
% Delete first "LAGS" rows to match the dimensions of X matrix
Y1 = Y1(pvar+1:Traw,:); % This is the final Y matrix used for the VAR

% Traw was the dimension of the initial data. T is the number of actual 
% time series observations of Y and X (we lose the p-lags)
T = Traw - pvar;

Y = Y1;
X = X1;
Z = Z1;


SIGMA = Q;

% First get Ordinary Least Squares (OLS) estimators
% A_OLS = inv(X'*X)*(X'*Y); % This is the matrix of regression coefficients
A_OLS = (X'*X)\(X'*Y); 
a_OLS = A_OLS(:);         % This is the vector of coefficients, i.e. it holds
                          % that a_OLS = vec(A_OLS)
SSE = (Y - X*A_OLS)'*(Y - X*A_OLS);
% SIGMA_OLS = SSE./(T-K);
SIGMA_OLS = SSE./(T);

%-----------------Prior hyperparameters for bvar model--------------------

n = K*M; % Total number of parameters (size of vector alpha)

% Define hyperparameters
    a_prior = 0*ones(n,1);   %<---- prior mean of alpha (parameter vector)
    V_prior = 100*eye(n);     %<---- prior variance of alpha
    
    if prior == 2
        a_prior = [1.35*eye(M); -0.4*eye(M)];
        a_prior = a_prior(:);
        V_prior = 0.001*eye(n);
    elseif prior == 3 % Minnesota
        a_prior = [zeros(1,M*constant); 1*eye(M); zeros(M,(M)*(pvar-1))'];
        a_prior = a_prior(:);

        % Hyperparameters on the Minnesota variance of alpha
        a_bar_1 = .2;
        a_bar_2 = .2;
        a_bar_3 = .2;
        decay = 2;

        % Now get residual variances of univariate p-lag autoregressions. Here
        % we just run the AR(p) model on each equation, ignoring the constant
        % and exogenous variables (if they have been specified for the original
        % VAR model)
        sigma_sq = zeros(M,1); % vector to store residual variances
        for i = 1:M
            % Create lags of dependent variable in i-th equation
            Ylag_i = mlag2(Yraw(:,i),pvar);
            Ylag_i = Ylag_i(pvar+1:Traw,:);
            % Dependent variable in i-th equation
            Y_i = Yraw(pvar+1:Traw,i);
            % OLS estimates of i-th equation
            alpha_i = inv2(Ylag_i'*Ylag_i)*(Ylag_i'*Y_i);
            sigma_sq(i,1) = (1./(T-pvar+1))*(Y_i - Ylag_i*alpha_i)'*(Y_i - Ylag_i*alpha_i);
        end

        % Now define prior hyperparameters.
        % Create an array of dimensions K x M, which will contain the K diagonal
        % elements of the covariance matrix, in each of the M equations.
        V_i = zeros(K,M);

        % index in each equation which are the own lags
        % index in each equation which are the own lags
        ind = zeros(M,pvar);
        for i=1:M
            ind(i,:) = constant+i:M:K;
        end
        for i = 1:M  % for each i-th equation
            for j = 1:K   % for each j-th RHS variable
                if constant==1 % if there is constant in the model use this code:
                    if j==1
                        V_i(j,i) = a_bar_3*sigma_sq(i,1); % variance on constant
                    elseif find(j==ind(i,:))>0
                        V_i(j,i) = a_bar_1./(ceil((j-constant)/M)^decay); % variance on own lags           
                    else
                        for kj=1:M
                            if find(j==ind(kj,:))>0
                                ll = kj;                   
                            end
                        end                 % variance on other lags   
                        V_i(j,i) = (a_bar_2*sigma_sq(i,1))./((ceil((j-constant)/M)^decay)*sigma_sq(ll,1));      
                    end
                    
                else  % if there is no constant in the model use this:
                    if find(j==ind(i,:))>0
                        V_i(j,i) = a_bar_1./(ceil(j/M)^decay); % variance on own lags
                    else
                        for kj=1:M
                            if find(j==ind(kj,:))>0
                                ll = kj;
                            end                        
                        end                 % variance on other lags  
                        V_i(j,i) = (a_bar_2*sigma_sq(i,1))./((ceil(j/M)^decay)*sigma_sq(ll,1));            
                    end
                end
            end
        end
    
        V_prior = diag(V_i(:));

    end
    
    
    % Hyperparameters on inv(SIGMA) ~ W(v_prior,inv(S_prior))
    v_prior = M;             %<---- prior Degrees of Freedom (DoF) of SIGMA
    S_prior = eye(M);            %<---- prior scale of SIGMA

%========================== Start Sampling ================================

%         %Check Stationarity (optional)
%         chck=-1;    % Draw rho but  ensure stationarity
%         while chck<0

if prior == 0
    
    %--------- Posterior hyperparameters of ALPHA and SIGMA with Diffuse Prior
    
    % Posterior of alpha|Data,SIGMA ~ Multi-T(kron(SSE,inv(X'X)),alpha_OLS,T-K)
    sparseZ = sparse(Z);
    invSigma = inv2(SIGMA);
    VARIANCE = kron(invSigma,speye(T));
    V_post = inv2(sparseZ'*VARIANCE*sparseZ);
    V_post = (V_post+V_post')./2;   
    a_post = V_post*(sparseZ'*VARIANCE*Y(:));
   
    
    alpha = a_post + chol(V_post)'*randn(n,1); % Draw of alpha    

    ALPHA = reshape(alpha,K,M); % Draw of ALPHA
    
    % Posterior of SIGMA|ALPHA,Data ~ iW(inv(S_post),v_post)
    v_post = T;
    S_post = (Y - X*ALPHA)'*(Y - X*ALPHA);
    SIGMA = inv2(wish(inv(S_post),v_post));% Draw SIGMA    
    
        if DiagOpt==1; 
        SIGMA=zeros(size(SIGMA));
        S_post = 0*S_post;
        EE = Y - X*A_OLS;
        for ii=1:M
            S_post(ii,ii)= EE(:,ii)'*EE(:,ii);
            SIGMA(ii,ii)=inv2(wish(inv(S_post(ii,ii)),v_post));% Draw SIGMA    
        end
        end
    
else
    
    % --------- Posterior hyperparameters of ALPHA and SIGMA with
    % --------- Independent Normal-Inverse Wishart Prior
    
    VARIANCE = kron(inv2(SIGMA),speye(T));
    V_post = inv2(inv2(V_prior) + Z'*VARIANCE*Z);
    a_post = V_post*(inv2(V_prior)*a_prior + Z'*VARIANCE*Y(:));
    alpha = a_post + chol(V_post)'*randn(n,1); % Draw of alpha
    
    ALPHA = reshape(alpha,K,M); % Draw of ALPHA
    
    % Posterior of SIGMA|ALPHA,Data ~ iW(inv(S_post),v_post)
    v_post = T;% + v_prior;
    S_post = (Y - X*ALPHA)'*(Y - X*ALPHA);%S_prior + (Y - X*ALPHA)'*(Y - X*ALPHA);
    SIGMA = inv2(wish(inv(S_post),v_post));% Draw SIGMA
    
        if DiagOpt==1; 
        SIGMA=zeros(size(SIGMA));
        S_post = 0*S_post;
        EE = Y - X*A_OLS;
        for ii=1:M
            S_post(ii,ii)=S_prior(ii,ii) + EE(:,ii)'*EE(:,ii);
            SIGMA(ii,ii)=inv2(wish(inv2(S_post(ii,ii)),v_post));% Draw SIGMA    
        end
        end
      
end

%     A = [ALPHA(constant+1:end,:)'; [eye(M*(pvar-1)) zeros(M*(pvar-1),M)]];
%         ee=max(abs(eig(A)));
%         if size(ee,1)==0 || ee<=1
%             chck=1;
%         end
%         end
        
resids = (Y - X*ALPHA);

if prior == -1
    
    ALPHA = reshape(a_OLS,K,M);
    SIGMA = SIGMA_OLS;
    resids = (Y - X*ALPHA);
   
end
        

end

